A Continuum Description of Vibrated Sand 
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The motion of a thin layer of granular material on a plate undergoing sinusoidal vibrations is 
considered. We develop equations of motion for the local thickness and the horizontal velocity of 
the layer. The driving comes from the violent impact of the grains on the plate. A linear stability 
theory reveals that the waves are excited non-resonantly, in contrast to the usual Faraday waves in 
liquids. Together with the experimentally observed continuum scaling, the model suggests a close 
connection between the neutral curve and the dispersion relation of the waves, which agrees quite 
well with experiments. For strong hysteresis we find localized oscillon solutions. 
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O ! I. INTRODUCTION 

■ Very little is known about the laws governing the macroscopic motion of granular materials, or for short, "sand". 
Most of our information comes from either experiment or microscopic molecular dynamics calculations. An accepted 
continuum description of sand, analogous to the equations of hydrodynamics, is missing. If such a description exists, 
it would help our understanding enormously, much in the same way hydrodynamics has dominated our understanding 

, of fluids. 

One major difficulty is that sand behaves very differently in different flow situations. If at rest or nearly so, sand 
behaves like a solid, and the packing of particles is very important Jj, while to reach a fluidized state the particles 
7— I ■ have to be shaken quite violently. In an interesting early paper, Haffg deals with grains in a nearly compact state. 

The opposite limit of low densities, where particle interactions are dominated by binary collisions, is known as "rapid 
granular flow" . Using methods of kinetic theory, this limit has received a great deal of attention ||-|| . It leads to 
complicated three-dimensional hydrodynamic equations with non-Newtonian transport coefficients that still need to 
be tested against experiment. 

Great interest has been stirred by recent experiments in which a thin layer of particles is placed on a plate undergoing 
sinusoidal vibrations [^|-^0). Above a certain vibration amplitude regular and irregular surface- wave patterns are 
excited subharmonically, i.e. the frequency of the waves is half that of the driving. The observed phenomena are 
very similar to Faraday waves excited in a periodically vibrated liquid layer (e.g. pl|). Typical experimental data are 
• • . the phase diagram of different wave patterns as a function of frequency and acceleration, and the wavelength of the 
patterns as a function of driving frequency at fixed acceleration. The data turn out to be independent of container 
size or shape, so the observed patterns represent an intrinsic property of the dynamics of vibrated sand. 

The experiments tri gge red a host of theoretical work in which a wide range of approaches has been used: molecular 
dynamics calculations fl4| , |l8||, simplified particle dynamics [^2[ , semi-continuum theories ^,^,0 , phenomenological 
Ginzburg-Landau models p5|, phenomenological coupled map models [ p6[ and order-parameter models p7|. Most 
notably, recent molecular dynamics simulations have reproduced the experimental results in great detail p8| . The 
comparisons between the continuum-type models and experiments have, however, not been very detailed. Their focus 
was mostly on reproducing the localized excitations of the layer ('oscillons') that have been observed experimentally 
p8| . Studies of more general order-parameter models J2{],|27| indicate, however, that such localized waves can also 
arise in non-granular systems and are therefore not the hallmark of these systems. In fact, similar excitations have 



been observed recently in Faraday experiments with shear-thinning clay suspensions 30 1 . 

The goal of this paper is to come up with a model for waves in vibrated sand that is based on physically accessible 
variables and is sufficiently realistic to allow a meaningful comparison and test with experiments. At the same time 
it should be simple enough to permit also analytical investigations, which often provide insights that are hard to 
obtain by numerical means. Based on the observation that the dispersion relation of the excited waves exhibits a 
continuum scaling in the small- frequency regime |3l|, 18| 20 1, we propose a continuum model. As indicated above, 



finding a continuum model based on the microscopic behavior of individual grains would be a formidable task. For 
part of the period, sand rests on the plate in a compacted state, while in the remainder of the period it is in free 
fall, leading to a fluidized state. In addition, little is known about boundary conditions near a solid wall fH32"[. 
Therefore, we adopt a purely phenomenological approach that is akin to a shallow-water theory, i.e. the vibrated 
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sand is described by its height and its horizontal velocity only. The proposed model differs, however, in significant 
aspects from a fluid-dynamical description. In particular, the unvibrated sand layer exhibits no oscillatory response. 
In contrast to the Faraday waves in liquids, the excited waves can therefore not be viewed as arising from the resonant 
driving of damped wave modes of the sand layer. 

The paper is organized as follows. In section [n] wc introduce a one-dimensional continuum model of vibrated 



sand. We discuss the driving mechanism and the physical significance of the parameters. In section |I| the model 
is linearized around a flat layer, to obtain analytical results for the onset of the waves and the dispersion relation 
between the excitation frequency and the wavenumber of the resulting surface waves. The results are compared with 
recent experiments and simulations by Bizon et al. ]18| ]. In the next section IV we turn to the nonlinear behavior of 
the model. First we investigate physical origins for the experimentally observed hysteresis in the onset of the waves. 
Tuning the parameters in the model to give strong hysteresis, we find localized oscillon solutions. 



II. THE MODEL 



Central to our model is the experimental observation that the position of the bottom of the layer of sand is closely 
modeled by the motion of a single, totally inelastic particle J|l|] . This is because all the energy is lost upon impact in 
the inelastic collisions between the grains. The force driving the patterns is thus proportional to the acceleration of 
the inelastic particle, minus the acceleration of gravity g. At each impact, this relative acceleration -y(f) is strongly 
peaked and its strength is related to the velocity of impact. 

We consider very thin layers and assume that they can be characterized by their thickness and mean horizontal 
velocity alone. For simplicity, we will only consider one-dimensional motion. We will not address the nonlinear 
pattern-selection problem that arises in two-dimensional patterns (e.g. stripes vs. squares). From mass and mo- 
mentum conservation considerations we arrive at equations quite similar to those of a fluid layer in the lubrication 
approximation, except for some crucial differences to be elaborated below. The equations are 

k+(vh)x = (Dxhs)., (la) 

vt + vv s = -iQE) 7 == -Bv+(D 2 v s ).. (lb) 
y 1 + h s 

To contrast all physical quantities from their dimensionless counterparts, they carry an overbar. Equation ( |Ta| ) 
comes from mass conservation, with an Edwards- Wilkinson diffusion term on the right, which describes the tumbling 
of grains atop one another ]33| ]. 

Equation (|lb|) expresses the momentum balance where v is the horizontal velocity integrated over the layer height 
h p4[ . It contains a driving term proportional to both the acceleration 7(f) relative to a freely falling reference frame 
as well as the slope hx, since particle motion gets started only if the surface is inclined. The denominator reflects 
the assumption that upon impact the sand grains are isotropically dispersed but only those scattered out of the layer 
contribute to the horizontal flux. Note that this ensures finite driving in the limit h s — > 00. Without the denominator 
we numerically found wave solutions which became progressively higher and more peaked, leading to an unphysical 
finite-time singularity. The second and third term on the right describe the internal friction due to vertical gradients 
in the velocity field and the viscous-like friction due to horizontal gradients, respectively. Since the vertical gradients 
are not resolved in this thin-layer approach they lead to a bulk damping term. 

Our model (Q) differs from the fluid problem in a number of ways. Through the assumption of random scattering 
of the grains upon impact the driving is nonlinear in h s . In contrast to liquids the sand layer lifts off the plate when 
the acceleration of the plate exceeds g. Thus, the acceleration term 7(F) vanishes over large parts of the cycle and 
is largest when the sand layer hits the plate. In the following we will assume that it consists mainly of a series of 
S- functions. There is no surface tension in the granular material. Instead an additional diffusion term appears in the 
equation for h. 

It should be emphasized that the effective friction and diffusion coefficients B, D\, and D2 contain implicitly the 
effect of the solid plate and the varying degree of fluidization present at different frequencies. It would be quite 
natural to assume that the friction between the sand and the plate arises only during the phases when the layer is 
very close to the plate. Since the velocity is not continuous through the <5-like impact of the layer on the plate a <5-like 
friction term is mathematically ill-defined. For simplicity we therefore smear out the friction over the whole period 
and expect that it can be modeled by an effective value of the coefficient B. We expect that due to the graininess of 
the material the dissipation due to the vertical gradients will increase with decreasing slope of the surface since the 
faster flowing grains near the surface are hindered by the slower grains underneath. This is not unlike the effect of an 
angle of repose. The coefficients D\ and D2 for particle diffusion and viscosity are expected to depend on the typical 
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velocity of the grains which is related to the impact velocity vq. They will increase with the typical velocity of the 
grains, since increasing velocity enhances the diffusive transport as well as the momentum exchange due to collisions 
[B5lp6[, Thus, in general we have 



Di,2 = D li2 (v , h, v), B = B(v ,h, h x , v). 



(2) 



Since the experiments are performed over a small range of accelerations the impact velocity is essentially given by 
the frequency. Thus, the dependence on vq implies an apparent frequency dependence of the coefficients. Note, a 
frequency dependence would also arise if an averaging procedure could be applied in which the basic equations are 
averaged over a period of the driving. It should be emphasized, however, that this is not the origin of the frequency 
dependence considered in this paper. 

We now make all quantities dimensionless using the filling height Jiq as a length scale and f = (ho/g) 1 ^ 2 as a time 
scale. The dimensionless driving j(t)/g depends then only on the dimensionless acceleration 



r = 



Auj 2 



fJ 



(3) 



of the plate, where A is the amplitude of the harmonic driving with frequency Co. The dimensionless dispersion 
relation has then the form q = q(u>,T,ho/d), where d is the grain diameter. The remarkable observation from the 
experimental data |31 18 pO] is that in the low-frequency regime the dispersion relation is independent of ho/d, i.e. of 



particle diameter, where ho/d varies between 3 and 14. 



III. LINEAR THEORY 



To understand the instability of the flat layer, we linearize in the dimensionless variables H and V , with h/ho = 
h = 1 + H and ^/(gho) 1 ^ 2 = v = V. The transport coefficients are evaluated at h = 1, v = 0. The linearized equations 
of motion can be transformed into a single wave equation for H: 



Htt = -BH t + (BDi 
(£>i + D 2 )H txx - 



-l{t))H xx + 
D!D 2 H X 

XXX • 



(4) 



To simplify the analytical calculation we assume that 7(f) consists only of a series of 5-shocks with period T = 2ir/u>, 



r(t) = vo Ht-jT), 



(5) 



and neglect the force on the layer during the subsequent time intervals in which the layer is in contact with the plate. 
Note that vq is the impact velocity, hence it can be written as vq — f(T)/u>. The curve f(T) is shown in Fig. |l|, see 
jnj. As T rises above 1, the layer begins to bounce and v increases with T. Between 3.3 < T < 3.7 the layer is locked 
into a state where it never rests on theplate and vq remains constant. Above T — 3.7 period doubling occurs and "/(t) 
can no longer be written in the form (|5|). In a straightforward generalization, two different periods with T\ +T 2 = 2T 
appear. 

8.0 




2.0 3.0 
Acceleration r 



FIG. 1. The dimensionless impact velocity of a completely inelastic ball on vibrating plate as function of the acceleration F. 
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In between shocks, the solution of (^) is H oc e at e lqx , where the dispersion relations 

a 1 = -D l q 2 , a 2 = -B-D 2 q 2 , (6) 
correspond to pure relaxation. At the point of impact, the (5-function imposes a jump condition 

H[ +) - Ht } = v H xx . (7) 
The height H itself is continuous. By making the general ansatz 

H n = (a„e CTl ( t ~" T ) + /3„e CT2(t -" T) ) e lqx 

and requiring a„+i = sa n and /3„+i = s[3 n for the eigenmode we obtain for the amplification s of the eigenmodc 

If 2 si - s 2 



s = p±Vp 2 - sis 2 , p=- [s 1 + s 2 -v q 2 , (8) 

with Si = e ' 17, and s 2 — e a2T . Since Sj < 1 (cf. eq.(^|)) and p < 1 the condition for instability, |s| > 1, can only 
be satisfied with real s < — 1. Thus, the only instability is subharmonic, i.e. the motion repeats itself only every 
other period of the driving, as observed experimentally. Analyzing d§|) in the limit q — > oo shows that catastrophic 
instabilities occur if either D\ or D 2 vanish. This can also be seen directly from ([!]): Since the driving with 7 is 
proportional to q 2 , only the combined damping D\D 2 q A keeps short-wavelength instabilities at bay. We emphasize 
that this mechanism for wavenumber selection is quite different from that at work in the liquid case. Faraday waves 
have the same dispersion relation as if there was no driving, the significance of the driving lies only in exciting the 
waves. In the case of sand, the medium itself does not support waves as seen from (|^); only the competition between 
driving at small wavenumbers and damping at large wavenumbers selects the wavenumber. 

Strictly speaking, (0) applies only during the free fall between shocks. During the periods in which the sand is in 
contact with the plate and experiences an acceleration it could exhibit oscillatory response if the damping is sufficiently 
weak. Specifically, in the absence of any forcing, i.e. for 7=1, the waves exhibit damped oscillatory behavior for 

B(D 2 - Di) < 1. (9) 

Experimentally, however, the layer spends a large fraction of the cycle in free flight already before the onset of waves. 
Therefore, even if (0) should be satisfied during the brief periods during which the layer is in contact with the plate 
the oscillatory response is not expected to be relevant for the excitation of waves. 

Next we compute the most unstable mode on the neutral curve s = 1. For supercritical transitions it gives the 
wavelength that is expected to appear as the acceleration is raised slowly above the critical value r c . Introducing the 
dimcnsionless combinations 

6=^, P = BT, Q 2 = D iq 2 T, (10) 



we find 



_ (1 + eQ 2 )(P +(S- 1)Q 2 ) 1 + e ~(^Q 2 ) 

1 Q2 I _ e _ (/3+ (5-l)Q2) ■ \ 1L ) 



To find the critical wavenumber, ( |TT|) has to be minimized with respect to Q, giving 

Q c = Qc(/3,5). (12) 
Plugging this back into ([ll]) leads to the critical impact velocity Vq(@, 5), and then to an expression of the form 

f(T c )=v^ ) uj = D 1 u;^(0,6). (13) 

In general, the minimum has to be found numerically. For two limiting cases, however, the dispersion relation can be 
given explicitly, 

q c =^^ 2 , for 00, (14) 

q c = a{D 1 ,D 2 ) w 1 / 2 , for/3^0. (15) 
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where a{D\,D2) is determined from an implicit equation. In both cases q c oc uj 1 / 2 for fixed Di. This is because the 
waves are damped by a diffusive mechanism. The same "viscous" scaling has also been found in other approaches 
^4j,^5|. Power laws different from 1/2 are due to some frequency dependence of the model parameters. 

We can now attempt to make a more quantitative comparison between experimental measurements fl8|| and our 
model. We will try to extract the dependence of the transport coefficients on the layer height and the frequency from 
the experimental data, and see whether this leads to a consistent picture. In Figure |2j and | we show the experimental 
measurement of the onset curve and the dispersion relation, respectively. In a double-logarithmic plot, the dispersion 
relation shows a sharp transition at oj = Lo tr , which lies between 3 and 4 for ho = 2.98mm. In the low-frequency 
regime u> < cutr the behavior is close to k ~ ui 1 - 3 , while a much weaker dependence of roughly k ~ ui ' 3 is seen at high 
frequencies. Note that the exponents of both power laws are much smaller than the exponent 2 reported in earlier 
experiments jllj, continuum theory J23J], and numerical simulations [|l4|,|l7|]. In view of the short scaling range one 
should, however, be very careful interpreting the data in terms of scaling laws. 

3.00 i i i i i i i i i i i 




Frequency co 

FIG. 2. Neutral curve T c (uj). Solid circles give the experimental results for increasing T for a mean layer thickness of 
2.98mm, taken from [uij. Theoretical fits for B = 0.08, Dj = 2.1/w, D 2 = 0.12/u (open squares), B = 0.14, Di = 0.93/w, 
D 2 = 0.94/w (plus), and B = 0.3, Di = 0.1/w, D 2 = 1.9/w (open circles). 
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Frequency co 

FIG. 3. Dispersion relation q c (u) corresponding to the neutral curves in Figj^. Experimental results are given by solid 
symbols (circles for ho — 2.98mm, triangles for ho = 1.49mm). 

According to |?l| the transition between the different power laws is associated with the frequency Ty/ho/2d, above 
which the velocity of the plate is no longer sufficient to let one particle hop across the other. Therefore, all particles 
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are locked into a fixed relative position and the motion is predominantly in the vertical direction. Our theory is 
therefore not expected to be applicable. Hence we will only be concerned with the low frequency regime U) < u>t r . 
From the data of |3l| for ho/d between 3 and 13 it is also seen that continuum scaling works much better in the 
low-frequency regime. 

Turning to the phase diagram, stripe patterns are observed at high frequencies while our frequency range of interest 
lo < Lu tr is associated with two-dimensional square patterns. This does not, however, affect the comparison with the 
linear properties of the one-dimensional model. In the low-frequency regime, the critical T decreases slightly with 



frequency, and hysteresis is found. Experiments at different layer heights |16 37]] reveal that there is only a small 
increase of r c with layer height. 

To compare theory with experiment, it is useful to return to dimensioned quantities in order to resurrect the 
dependence on the mean layer height Tiq. We find 

f(T c ) = ^^-u>mS), (16) 

t = Q - M6) : *• (17) 
27rL>i(w,/i ) 

Assuming that both diffusion coefficients scale the same way in ho and u>, we try the ansatz 

/> 

Dx,2 =£>i,2-r£- (18) 

This renders 5 independent of u> and h$. Since experimentally T c increases only slightly with ho at fixed u>, we 
conclude from ([l6]) that /x = 1 + e, where < e <C 1 is used to indicate the weak increase in T c with layer height 
p7| . This is not to say that the experimental onset fits such a power-law, but rather just to indicate how a change 
in the dependence of the onset on ho affects the dispersion relation between q and to within the present framework. 
In all results below we use e = 0. Inserting ( |l8|) into the dispersion relation ( |l7j ) one obtains for the dimensionless 
wavenumber q as a function of lo 

f= & > 5 l h^^. (19) 

Experimentally, q(co) collapses onto a single curve independent of ho- Thus v = 1 — 2e and we obtain 

h 1+t 

^1,2 = ^1,2^7, ( 20 ) 

where Z?i.2 has the dimension of an acceleration (for e = 0). This makes the dimension of the material parameters of 
the sand the same as the control parameter of the experiment, which is the physical origin of the observed collapse 
of the dispersion relation. Note that (for e = 0) D12 are proportional to the impact velocity vq oc lo^ 1 for given 
acceleration. Thus, the final result for the dispersion relation is 

^ = «MC^). (21) 

It remains to determine 0, S and D\ from a comparison with the experiments. 

From (|l^) and (|2C]) it follows that for B = the onset acceleration T c is independent of frequency (for e = 0). If 
B is taken to be finite, lower frequencies will be damped most, due to the term Bv in (|TJ). Thus, in accordance with 
experimental findings |l2|JT^ , p^ |, the critical acceleration r( c ) decreases slightly with frequency. 

We concentrate on the transition from a flat layer to waves in the low- frequency regime lo < Lo tr . The transition 
occurs for T between 2.5 and 3 at the lowest frequencies, depending on the type of particle and on the layer thickness. 
For definiteness, we consider the phase diagram jL8) for a mean layer thickness of ho — 2.98mm, shown in Fig.|[ The 
critical T is seen to rise from r c = 2.3 at u> = 4 to T c — 2.7 at lo = 1. Keeping the ratio S — D2/-D1 constant, one can 
adjust D2 to give the correct value for lo = 1. For B — 0, this value would remain constant for all lo. By choosing B 
to have F c agree with the experimental finding at lo = 4, the onset curve is reasonably well reproduced. In Figure |2| 
we include the fits for S = 0.057, S = 1.01, and S = 19, for which the agreement is comparable. 

Next we turn to the dispersion relation of wavenumber versus frequency as given by Bizon et al. Jig ], shown in 
Figure ||. The only parameter remaining to be fixed is the ratio S between damping constants. We adjust it to obtain 
an optimum fit of the theoretical dispersion relation to the experimental data. While the slope of the experimental 
data is still slightly higher than predicted from theory, for S — 0.057 we find a reasonable agreement. 
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IV. NONLINEAR PROPERTIES AND OSCILLONS 



To investigate the nonlinear behavior of the model, the coefficients Di are taken to depend on the local layer 
height h according to ((2^). We find that the bifurcation from the flat state to the standing waves is supercritical if 
the transport coefficients are taken to be independent of h x and v (as in the similar model proposed in p3|]), while 
experimentally the transition is subcritical |ll| ,jf2| . We consider two physically plausible reasons for this discrepancy. 

Statically, sand has a finite angle of repose. Motion only starts if the slope of a hill of sand exceeds a certain critical 
value k. One may expect that in a dynamic state this leads to enhanced friction of the flowing upper layer when the 
slope of the surface is small. We model this by taking 



B = B Q [1 + Bte 



-(K/k)'- 



(22) 



where k sets the characteristic slope below which the granular character of the material becomes noticeable. 

Second, the layer will have a much higher viscosity when it is near its compact state. An accurate description will 
have to contain a parameter which measures how far the layer is from this state. As the "fluidization" increases, 
the viscosity is expected to decrease. The simplest assumption is that the fluidization depends directly on the local 
velocity v(z, t), and thus we model this effect by 



D 2 =D^ ( W ) 



f + rye 



-(v/v f ) 



(23) 



Note that the fluidization and thus the local velocity is distinct from the typical velocity of the vibrating plate. The 
latter sets the typical collision time between particles and thus leads to a rise in viscosity as it increases. It is obvious 
that (|22] ) and (|23j ) can lead to a subcritical transition to waves. Greater accelerations are needed to set the layer into 
motion for zero initial slope and zero velocity. On the other hand, once waves have started to appear, the motion will 
persist to lower values of T. Fig. || shows a typical hysteresis with only the fluidization ( |23| ) taken into account. In 
the numerical simulations j(t) is taken to contain not only the <5-shocks but also the smoothly varying acceleration 
while the layer is in contact with the plate. 
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FIG. 4. Hysteretic nature of the transition in the presence of a velocity-dependent viscosity (|23|) with r\ — 0.5 and Vf = 0.2. 
The other parameters are u = 1.62105 Di = D 2 = 0.37, and B = 0.45. 

The subcritical properties of the model are believed to play a crucial role in the appearance of localized subharmonic 
excitations called oscillons. They arise from extended square patterns when the driving is reduced below the stability 
limit of the latter. Alternatively, they can be excited by a localized "seed" j28|. In one period, an oscillon consists of 
an axisymmetric jet of a height several times its diameter shooting out of an almost flat layer. In the next period, 
the oscillon forms a shallow circular trough, surrounded by a small mound. 

The appearance of oscillons requires that the flat layer be linearly stable, while in the region covered by the oscillon 
waves can be sustained, since the grains are fluidized only in the center. In the "trough" -phase, the mound has a 
very shallow slope facing out, so the outward flux of material is strongly damped and almost all the material is sent 
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radially inward. The small amount of sand that is transported outward comes back during the following "peak" -phase. 
Because of mass conservation, it is evident that the formation of oscillons is greatly facilitated by their radial geometry. 
The greater the radius r away from the center, the smaller the height h(r) that corresponds to the mass of the peak 
during the "peak"-phase. Conversely, ingoing sand is focussed into a sharp jet. Therefore, within the one-dimensional 
version of our model oscillons are expected to arise only over a smaller range of parameters than in two dimensions 
f27f . In fact, they could only be found for greatly exaggerated subcritical behavior, i.e. large B\ and r\. This may also 
be the reason why so far stable oscillons have not been found in experiments on vibrated sand between two narrowly 
spaced plates, which effectively have only one horizontal dimension Jig] , In these experiments localized structures 
appear only intermittently as localized bursts, but not as time-periodic structures with steady amplitude. Our main 
point will therefore be that structurally our model allows oscillon solutions. For a more quantitative description an 
axisymmetric version of our code has to be considered. 
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FIG. 5. Oscillon solution for T = 2.5, uj = 1.6, Di = 0.376/i, D\ " = 0.08/i, Bo = 0.1, Bi = 19, k = 0.2, r) = 6, «/ = 0.4. 

Figure [| shows a numerically obtained one-dimensional oscillon in its two phases. The central jet is considerably 
less sharp than the experimental one, as expected from the above argument. The parameters are chosen so as to damp 
the motion outside the oscillon very rapidly (r? = 6) and to make sure that as the layer hits the plate in the trough 
phase, almost all material is transported inward (k = 0.2, B\ = 19). The temporal evolution of the surface height is 
shown in Fig.^ as a space-time diagram. It illustrates the transport away and towards the center of the oscillon. As 
initial conditions, we chose a localized velocity profile directed inward towards a point, to produce an initial central 
peak. We find that the layer thickness averaged over two periods is smaller inside the oscillon than outside. Thus, 
the oscillon tends to push out material. In preparing an initial condition, we accounted for this by reducing the layer 
thickness in the center. Once the solution was close to stationary, we checked for exponential convergence towards 
an oscillon solution. Furthermore, the parameters could be varied to within a few percent without destabilizing the 
oscillon. 
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FIG. 6. Space-time diagram of the evolution of the oscillon shown in Fig.|5| during two periods of the driving. The shock-like 
impact occurs at the times marked with narrower line spacing. 

V. DISCUSSION 

In the present paper we have presented a simple hydrodynamic description of surface waves on vertically vibrated 
granular media. We have eliminated all vertical degrees of freedom and describe the horizontal motion by three 
effective transport coefficients. Our model predicts a new type of instability based on the interplay between shocks 
imparted by impacts on the plate and diffusion both in the layer height and the momentum of the layer. This 
distinguishes the model from earlier approaches |Tl|Jl6| , p5t , which in analogy to Faraday waves in liquids are based on 
the resonant excitation of damped waves that exist even in the absence of periodic driving. 

To obtain quantitative agreement with experiment, the dependence of the transport coefficients B, D\, and D2 
on experimental control parameters has to be taken into account. Remarkably, the possibilities for this dependence 
become severely restricted by the weak dependence of the neutral curve on the frequency and the collapse of the 
experimental data for the dispersion relation in units of the gravitational acceleration g and the mean layer height 
ho. As a consequence, the neutral curve and the fact of the data collapse alone imply a certain dependence of the 
wavenumber on the frequency, which turns out to be in reasonable agreement with the experiment. This leads to an 
interdependence of these seemingly independent observations. Conversely, the power laws are predicted to change if 
other parameter dependencies of B, Di, and D 2 apply, in which case data collapse may also no longer occur. Thus, the 
data collapse indicates more than continuum scaling alone, i.e. more than the independence of the grain diameter d; 
within our model the collapse is due to specifi c dy namical properties of the vibrated sand as reflected in the transport 
coefficients. Indeed, there are indications [jl6|,|l8[ of a high-frequency regime in which the power law is different and 
where also the data collapse seems to be in question pl[ . 

To obtain the experimentally observed strong sub-criticality of the transition to waves we have invoked a critical 
slope k in the friction term (p^). This effect would be connected with the granularity of the material and the 
parameters are expected to depend on d/hp, implying that the continuum scaling should not hold in the nonlinear 
regime. Unfortunately, the fluidization ( ]25| ) also comes into play when determining the hysteresis, so there is no simple 
estimate of the dependence of the hysteresis Ar on the particle diameter. Still, measurements of Ar as a function of 
both d and ho would be revealing. 

It is expected that the agreement with experimental data could be further improved if additional aspects of the 
system are included in the model. For instance, we have neglected effects of a layer dilation, which will depend on 
frequency; the dilation will smooth out the <5-driving, making the driving less effective, and lowering Vq- 



It would be interesting if the dependence of -B, D%, and D% on the frequency or other experimental variables could 
be determined directly in experiments. Whether the measurement of the frequency dependence of the viscosity by 
dragging a sphere through a vibrated layer of sand p3] is directly transferable to our model is not clear. In addition, 
it is assumed that only a certain part of the layer takes part in the horizontal motion. In [^| this was treated 
by introducing a "penetration depth" , which is an unknown function of parameters. It was crucial for obtaining a 
transition at io tr between a high and a low frequency regime pj] . 

An obvious question is whether our model can be tested against other experimental observations without signifi- 
cant further complications and with the same parameter values. First, the success |18| of numerical simulations in 
reproducing experimental data opens the possibility of a detailed comparison of the wave forms. Nonlinear effects 
are very strong and the waves are far from sinusoidal, the crests being substantially more peaked than the troughs, 
which qualitatively is also seen in our model. A next step would be an exploration of the phase diagram at higher 
accelerations, at least in the low frequency regime. That would necessitate a two dimensional approach to address 
the relevant pattern-selection properties. Such a two-dimensional version would also allow for a quantitative study of 
oscillons, both in their shape and their occurrence in parameter space. 

In the model presented here it is found that localized waves expel material into the quiescent regions. Transport of 
this type was considered essential in the Ginzburg-Landau-type model introduced in p5fl , in which the layer height 
couples back to the damping of the waves. Such a feedback would occur in the model discussed here for e > in (|2C|). 
Whether such a transport actually occurs in the experiment is of great interest. 

A significant step beyond the models presented so far would be an investigation of the high-frequency regime where 
the dispersion relation very significantly changes its power law. Since this seems to be the regime where grains 
are geometrically obstructing each other, one expects a significant dependence on grain size and a breakdown of 
continuum scaling. Within our model, this could possibly be modeled by a different parameter dependence of the 
transport coefficients. On the other hand, it is conceivable that the strong spatial correlations introduced by the 
particles, which are locked into a fixed relative position and predominantly move in a vertical direction, calls for a 
different, perhaps nonlocal hydrodynamic description. 
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